A stable Algebraic Spin Liquid in a Hubbard model 
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We show the existence of a stable Algebraic Spin Liquid (ASL) phase in a Hubbard model defined 
on a honeycomb lattice with spin-dependent hopping that breaks time-reversal symmetry. The 
effective spin model is the Kitaev model for large on-site repulsion. The gaplessness of the emergent 
Majorana fermions is protected by the time reversal (TR) invariance of this model. We prove that 
the effective spin model is TR invariant in the entire Mott phase thus ensuring the stability of 
the ASL. The model can be physically realized in cold atom systems and we propose experimental 
signals of the ASL. 



The concept of a spin liquid as a Mott phase with- 
out any local order was put forward by Anderson 
PP. Rs relevance to the physics of high temperature 
superconductors [H [3] led to the development of a gauge 
theory of spin liquids [3l [4] , analogous to quantum elec- 
trodynamics(QED). The spinons are the counterpart of 
electrons in QED and the visons, another emergent exci- 
tation, are the counterpart of the photon. Attempts at 
understanding the emergence of fermionic quasi-particlcs 
in spin systems in analogy with the anyonic quasi- 
particles in fractional quantum Hall systems has led to 
a general theory of quantum/topological order in spin 
liquids 5 . Experimental evidence for a spin liquid ground 
state has been seen, for instance, in the organic material 
K-(BEDT-TTF) 2 Cu 2 (CN)3[g. 

Frustration in magnetic interactions and quantum 
fluctuations tend to prevent magnetic ordering. Thus 
ASLs have primarily been studied in frustrated spin- 
1/2 Heisenberg anti-ferromagnets. [7]-[9] . The ASL shows 
power law decay not only for spin correlations but many 
other local order parameters as well. Hence it is intrinsi- 
cally susceptible to one of them ordering and inducing a 
spinon gap. Thus any realization of this phase must be 
accompanied by a mechanism of ensuring its stability. 

Kitaev[TU] constructed an exactly solvable anisotropic 
spin-1/2 model on a honeycomb lattice that exhibits the 
important properties of an ASL. It can be expressed as a 
model of two gapless Majorana-Dirac fermions (spinons) 
interacting with Z2 gauge fields (visons). A remarkable 
feature of the model is that the magnetic flux associated 
with every plaquette is conserved, as a result of which the 
visons are static. Consequently, while multi-spin opera- 
tors that conserve flux have algebraic correlations, those 
which do not, including the single spin operators, are ex- 
tremely short ranged [TT]. However, it has been shown 
that, for a certain class of perturbations that break these 
conservation laws, the spin-spin correlations become al- 
gebraic as well 8 , 9 . ASLs are thus realized in such per- 
turbed models. The stability of the ASL in the Kitaev 
model is due to time reversal (TR) symmetry: the two 
Majorana-Dirac fermions combine to form a single Dirac 



fermion, with an energy spectrum that cannot have a 
gap without breaking TR symmetry. An exactly solv- 
able spin-3/2 model with algebraic spin correlations has 
also been constructed [H]. 

The ASL has not yet been realized in a model of inter- 
acting fermions, though attempts have been made spec- 
ulating the possible means of doing so fTSrfT?] . However, 
it was shown [T5] that a short-ranged spin liquid (SRSL) 
emerges between the semi-metal and the Neel phases in a 
Hubbard model defined on a honeycomb lattice. A recent 
work claims otherwise [19]. 

The above discussion suggests that a Hubbard model, 
which is described effectively by the Kitaev honeycomb 
spin model in the large U limit, is a good candidate for 
realizing an ASL. Such a model was proposed by Duan et 
al. as a way of realizing the Kitaev model [10] in cold atom 
systems [20]. This model, which we henceforth call the 
Kitaev-Hubbard model, has anisotropic spin-dependent 
hopping which leads to the high degree of frustration in 
the effective spin model. The Hamiltonian of this model 
is 
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where Cj CT annihilates a fermion of spin projection a =t, i 
at site i (the spin index is implicit in the first term), a a 
(a = x, y, z) are the Pauli matrices, n a = c\c a is the 
number of fermions of spin a at site i, and hij) a denotes 
the nearest-neighbor pairs in the three hopping directions 



of the lattice (see Fig. SI I. 

At t' = 0, the model reduces to the simple spin- and 
TR-invariant, nearest-neighbor Hubbard model[18l [L9] . 
The term proportional to t' is a spin-dependent hopping 
term and breaks TR symmetry, SU (2) spin symmetry 
and the three-fold spatial rotation symmetry of the t' = 
model. It is however invariant under a spatial rotation 
of 2n/3 combined with a spin rotation of 27r/3 about 
the (111) spin axis. At t' = t, the one-body part of 
the Hamiltonian is a combination of the projection op- 
erators |(1 + <j a ). Thus, only those electrons that are 
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FIG. SI: (color online) The honeycomb lattice with the two 
sub-lattices marked by white and black dots. The six-site clus- 
ter used in this work is shown as the shaded area. The en la- 
bel the different spin-dependent hopping directions (blue solid 
lines), whereas the inter-cluster bonds are shown as dashed 
lines. 
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FIG. S2: (color online) The phase diagram of the Kitaev- 
Hubbard model at half-filling, showing the phases. The tran- 
sition from the AFI to the ASL phase is discontinuous. The 
red squares correspond to the parameter values at which the 
spectral graphs have been plotted in Fig. [S3] 



spin-polarized in the a th direction can hop along the a 
bonds. At this value of t', the effective low-energy spin 
model, at half-filling and large U, is the Kitaev honey- 
comb model. [201 [21] 

At U = and t' = 1 the single particle spectrum of this 
model shows four distinct bands, each of which has a non 
zero Chern number v\12\. The top and the bottom bands 
have v = 1, while the two middle bands, with v = — 1, 
are connected at the Dirac points. At t' = 1, the top 
two as well as the bottom two bands are gapped. As t' 
is decreased, this gap shrinks and finally closes at t' — 
0.717. The existence and locations of Dirac points can be 
experimentally measured in optical lattice systems 23 . 

As mentioned earlier, at t' = t, U — > oo, the model is 
analytically tractable and as we show later, is an ASL. We 
have computed the extent of the ASL phase using Clus- 
ter Perturbation TheorypQ (CPT) and the Variational 
Cluster Approximation (VCA).[3] These numerical tech- 
niques allow us to map the charge-gap of the model on to 



FIG. S3: Spectral functions of the Kitaev-Hubbard model, 
computed using CPT, as a function of energy (cj, y-axis) and 
momentum (k, x-axis) for the four sets of parameter values 
marked by the squares in Fig. |S2| red indicating maximum 
value and blue indicating minimum values. The spectrum 
is gapless only for the SM. F, M and K represent the high 
symmetry points of the Brillouin zone. 



the t' — U plane, and to calculate the extent of the Neel 
phase. VCA also allows us to find out whether or not the 
transitions out of the Neel phase are continuous. How- 
ever, the same cannot be done for the Mott transitions 
to the spin-liquid phases for which a cluster dynamical 
mean field technique would be required. 

Setting t = 1, we study the phase diagram of the 
Kitaev-Hubbard model from t' = 1 to t' = 0.5. Our 
results are summarized in Fig. |S2| At low U there is a 
TR-breaking semi- metallic phase (SM), characterized by 
charged, spin-1/2 fermionic quasi-particles. This phase 
exists in the region 1 > t' > 0.5 and U < 2.4. When 
U « 2.4 and if = 1, there is a Mott transition from the 
SM to the ASL phase, which then extends for large val- 
ues of U at this t'. Between U > 1.5 and U < 2.4 and 
with steadily decreasing t' , the system starts off in the 
SM phase, then makes a transition into the ASL phase 
and finally re-enters the SM phase, wherein it continues 
to live till t' = 0.5. At values of U greater than 2.4, de- 
creasing t' destabilizes the ASL phase and brings about a 
transition to the anti-ferromagnetic (Neel) (AFI) phase. 
The width of the ASL phase in the t' variable is maxi- 
mum at U ~ 2.4, extending from t' = 1 to if w 0.7. 

The ASL phase is bounded by the AFI and SM phases 
and is hence not connected with the possible short-ranged 
spin liquid at t' = OpSQI]. The SM and AFI phases do 
not have quasi-particles with fractional quantum num- 
bers or statistics. Thus the ASL is topologically distinct 
from the SM and the AFI and we expect the transitions 
between them to be discontinuous. 
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Using CPT, we obtain spectral functions correspond- 
ing to each of the three phases and plot them in Fig. S3 
Our numerical computations show that the transition be- 
tween the AFI and ASL phases is discontinuous (first- 
order) as expected. It is intriguing that the single parti- 
cle bands of opposite Chcrn numbers remain gapped in 
the same range of values of t' as the existence of the ASL. 
This seems to indicate that geometric phase effects may 
play an important role in this model. 

We study the ASL spin-spin correlations by deriving 
an effective, large-/7 spin Hamiltonian of the Kitaev- 
Hubbard model at half filling. To leading order in 1/U, 
this effective Hamiltonian is 
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This is a combination of the Heisenberg and Kitaev mod- 



els. The effective Hamiltonian ( S2 1 is TR symmetric, and 



this is thus responsible for the gaplessness of the Majo- 
rana spinons. Using charge-conjugation symmetry, we 
have analytically proved that the effective spin model is 
TR symmetric to all orders in 1/U (see Supplementary 
Information). Therefore, the spinons remain gapless in 
the Mott phase. As mentioned earlier, in the Kitaev 
model, conservation laws cut off the spin-spin correla- 
tions beyond nearest neighbors [TT]. Therefore, in order 
to analyze these correlations, we need to consider higher 
order terms. 

We have calculated the next order term in the effective 
spin Hamiltonian and obtained 
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where Lij) a denote nearest neighbors in the a direction 
and LiLdj)) a p denote next-nearest neighbors reached by 
first moving in the a direction and then in the /3 direction. 
The spin-spin correlation is, 



g(r,t) = (T (SS(<)5^(0) 



(S4) 



where r is a site of the Bravais lattice, I is the sublattice 
index. We compute this perturbatively in the fermionic 
representation of the spins, |7] in which the Majorana 
fermion operators Cj and bf are defined as follows: 



a? = iabf, 
{b<*M} = 25 ct f j 8 ij , 



{ci,b?} = 



(S5) 



The physical subspace is defined by the constraint 

WMphys = l^>ph ys (S6) 

In terms of these Majorana fermions, the leading order 
Hamiltonian is, 
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where J 
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This Hamiltonian describes 



V u 3 

Majorana fermions (cj), which we refer to as spinons, 
propagating in the background of static Zi gauge fields 
( u (ij)a = ibfb'j). The spin-spin correlation functions 
therefore factorize into propagators of the Ci operators. 
Since the spin operators create two units of flux on ad- 
joining plaquettes, the Majorana fermion propagators are 
in the background of an even number of fluxes at a few 
points. 

We are interested in the asymptotic form of the leading 
order correction, g( 2 \r,t) : when |r|, t — > oo. Tikhonov 
et. al. [5] have shown that in this limit, the propagators 
are the same as those in the flux-free background. Their 
result can be physically understood by noting that the 
particles hopping far way from the flux pairs will not 
pick up any phases from them. Thus we can expect the 
long wavelength modes to be insensitive to a few localized 
flux pairs. 

To compute the asymptotic form or the propagators, 
we can derive the continuum theory of the low-energy 
modes in the flux-free background. The Hamiltonian 
(S47), when ibfbj = 1, reduces to nearest-neighbor 
hopping on a honeycomb lattice just as in graphene. 
Graphene has low-energy Dirac quasi-particles about two 
points, K and K', in the Brillouin zone. However, since 
the c fermions are Majorana fermions, the excitations ex- 
ist only over half the Brillouin zone. Thus the low-energy 
modes constitute a single Dirac quasi-particle. The con- 
tinuum theory is derived by introducing slowly varying 
fields ipi (r) such that 
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ip(r, t) satisfies the gapless Dirac equation. Any mass 
term for a single Dirac fermion breaks TR invariance. 
Thus the gaplessness is protected by the TR invariance 
of the effective spin model. 

We have calculated the dynamical spin-spin correla- 
tions in the long-wavelength limit following the methods 
of Tikhonov et aZ.,[Sj (see Supplementary Information): 

(s z rl {t)s z M) 

= (-0.56 cos(2K • r)jf + I.I37172 + 0.887^) det G 
+ 0.28^ 2 ti(r x GT x G) + (0.56 7l 2 + 0.16^ 2 MGt x G) 



(O.287? + O.O27! + 0.16 7 i7 2 )(trG) 2 
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where 71 = 72 
(27r/3,27r/3^) and 
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Here t is the time, t = (t 1 
trices. 



and K is the Dirac point 



(S10) 



47T ( r 2_ J2 t 2); 
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Using equation (S51 ) and (S10), we find that the long- 
wavelength correlation function falls off as 1/r 4 . This 
algebraic decay has also been seen for single-spin per- 
turbations studied in Ref. |S1 Although the prefactor is 
extremely small (~ 1/C/ 6 ), this is the leading behavior at 
long distances. Therefore the effect of the perturbation 
cannot be neglected for any value of U, however large. 
Indeed, we can expect the strength of these correlations 
to grow as U decreases. This proves the existence of the 
ASL in the Kitaev-Hubbard model. 

Thus, at large U, the leading order contribution to 
the spin susceptibility is independent of U as in the Ki- 
taev model, whereas the next order contribution goes as 
(t/U) 6 . The U dependence of the spin susceptibility will 
hence be of the form \ = a + b(t/U) 6 , where a and b 
are constants independent of U. Experimental methods 
for measuring the spin susceptibility in cold atom sys- 
tems have recently been developed. [27] Thus susceptibil- 
ity measurements as a function of U can provide evidence 
for the existence of the ASL in this model. 

The stability of the ASL in this model comes from the 
preservation of TR symmetry We have investigated the 
possibility of spontaneous breaking of TR symmetry and 
the consequent emergence of a chiral spin liquid (CSL) 
with a spinon gap (see Supplementary Information). In 
the fermionized version [TU] of the effective spin model, 
where H e g = + we use a mean-field theory 

in which the vison and the spinon sectors are decoupled. 
This mean-field theory is exact for the Kitaev model, 
which is obtained by putting t 1 = 1 in . We find 
that the CSL solutions occur only for U < 1.6 and for 
0.5 < t' < l,and thus are not seen in the Mott regime 
U > 2.4. In Fig. S7 we plot the spinon gap as a function 
of U for t' = 1. 

In conclusion, we have shown that the Kitaev-Hubbard 
model, shows a Mott transition from a semi-metallic 
phase, to an algebraic spin-liquid phase. The former 
breaks time-reversal symmetry whereas the latter pre- 
serves it. The ASL is stabilized by TR symmetry. We 
have proved the TR invariance in the Mott phase (to all 
orders in t/U), using charge-conjugation symmetry. At 
intermediate U the ASL phase occurs for a wide range 
of t' which narrows down as U is increased. Concrete 
schemes to realize this model have been proposed [2TJ1 12"5] . 
and experimental methods to probe the semi-metal at 
low U 23 and the ASL at large C/j27j exist. This demon- 
stration of the existence of the ASL might help better 
understand the physics of the pseudo-gap phase of the 
underdoped high termpertaure superconductors [51 [51 12"5] . 
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FIG. S4: (color online) Gap of the spinon spectrum as a func- 
tion of U for t' = 1 
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Supplementary Information 
Cluster Perturbation theory and the VCA 

Cluster Perturbation Theory (CPT) is an approxima- 
tion scheme for the one-electron Green function G(cj) 
within Hubbard-like models. [IH3] It proceeds by dividing 
the infinite lattice 7 into a super-lattice T of identical 
clusters of L sites each (Fig. 1 of the main text illustrates 
the cluster used in this work). The lattice Hamiltonian 
H is written as H = H c + Ht, where H c is the cluster 
Hamiltonian, obtained by severing the hopping terms be- 
tween different clusters, which are put into Ht- Let T 
be the matrix of inter-cluster hopping terms and G c (uj) 
the exact Green function of the cluster. Because of the 
periodicity of the super-lattice, T can be expressed as a 
function of the reduced wave-vector k and as a matrix in 
site indices within the cluster: T a b(k). Likewise, G c is a 
matrix in cluster site indices only, since all clusters are 
identical: G c ab {uj). Thus, hopping matrices and Green 
functions in what follows will be k-dependent matrices 
of order L, the number of sites within each cluster. The 
fundamental result of CPT for the system's one-electron 
Green function is 

G- 1 (k,a;) = G c - 1 (w)-T(k) (Sll) 

In practice G c (uj) is calculated numerically by the Lanc- 
zos method and the cluster must be small enough for this 
to be possible. Because the lattice tiling breaks the orig- 
inal translation invariance of the lattice, a prescription 
is needed to restore the translation invariance of the re- 
sulting Green function. The CPT prescription for this 
periodization is 

G(k, u ) = i £ e-W-^-^G ab (k, u ) (S12) 

a, b 

where now k belongs to the Brillouin zone of the original 
lattice and the sum is carried over cluster sites. This 
formula is exact in both the strong (t — > 0) and the weak 
(U —> 0) coupling limits. 

Once the approximate interacting Green function 
can be calculated, the spectral function A(k, ui) = 
— 2ImG(k,w) follows. From there the density of states 
N(w) can be calculated by numerically integrating 
A(k, uj) over wave- vectors and the existence of a spectral 
gap can be assessed. Numerically, the density of states 
is always evaluated at a complex frequency with a small 
imaginary part r\ that broadens the spectral peaks. By 
applying a few values of r\ and extrapolating to 7/ — s- 0, 
one can detect the presence (or not) of a spectral gap at 
the Fermi level. This allows us to distinguish between a 
metal and a Mott insulator. 

If we know at what wave-vector the gap first opens 
up, as is the case at t' = (The Dirac points), then we 
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can estimate the gap more reliably, without the need to 
extrapolate to rj —¥ 0, by simply looking up the Lehmann 
representation of the CPT Green function, which can be 
calculated when the cluster Green function is computed 
using the band Lanczos method. 

The Variational Cluster Approximation (VCA) is an 
extension of CPT in which parameters of the cluster 
Hamiltonian H c may be treated variationally, according 
to Potthoff's Self-Energy Functional Theory (SFT).[1E>] 
In particular, it allows the emergence of spontaneously 
broken symmetries and provides an approximate value 
for the system's grand potential fl. In the case at hand, 



a single variational parameter is used: the strength M c 
of a staggered magnetization field that is added to the 
cluster Hamiltonian: 



H M = M^2m a clc 
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where the symbol m a is +1 for spin-up orbitals on the 
A sublattice and spin-down orbitals on the B sublattice, 
and —1 otherwise. 

Technically, VCA proceeds by minimizing the following 
quantity: 
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where VL C (M C ) is the grand potential of the cluster alone 
(obtained in the exact diagonalization process). The in- 
tegral over frequencies is carried over the positive imag- 
inary axis. At the optimal value M*, f2(M*) is the best 
estimate of the system's grand potential. At this value 
of M c , the order parameter M is calculated: 
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where G aa are the diagonal elements of the CPT Green 



function (Sill 



VCA provides estimates of order parameters, much like 
mean-field theory, but is quite superior to it because the 
Hamiltonian remains fully interacting (no factorization 
of the interaction) and spatial correlations are treated 
exactly within the cluster. 

In this work, the VCA was used to find the phase 
boundary of the antiferromagnetic phase. On the other 
hand, the transition between the spin liquid and semi- 
metal phases was found by monitoring the closure of the 
gap via CPT only. 

Proof of time reversal symmetry in the Mott Phase 



Let P m represent the projector on the Hilbert subspace 
containing m doubly occupied sites. Thus P^ = P m and 
V°° P = 1 

We write the Hamiltonian as 



oo 

T s = Pa+rnH-IcPm 
m=0 



(S19) 
(S20) 



s = -1,0,1, P_x = 0, T_ s = Tj and T s increases the 
doubly occupied sites by s. Thus we have 



[H ,T S ] =sUT s 



(S21) 



The operators that commute with T-Lq are called block 
diagonal (BD) and the others off block diagonal (OBD). 



In this section we will outline the derivation of the ef- 
fective spin Hamiltonian and prove that it is time reversal 
(TR) invariant in the Mott phase. 

Notation 



The Hamiltonian is taken to be: 
H = Ho+Hk 

Ho = nt^nii 



U K = £ 
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(S18) 



Canonical Perturbation Theory 



We follow the approach of Chcrnyshev et a/.[6 and 
block diagonalize the Hamiltonian order by order in t/U. 
The (k + l) th order Hamiltonian is written as 
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is chosen so as to eliminate the OBD terms of order 
(£) that survive in the Hamiltonian after performing 
the canonical transformation at order k — 1. By con- 
struction, £W does not contain terms that preserve the 
number of doubly occupied sites. 
S( k > has to satisfy the equation 
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where %ob D is the OBD part of . It is easy to show 



that 



(S24) 



The effective spin Hamiltonian at half-filling, H^ k \ is 
obtained by projecting H^- ' onto the singly occupied sub- 
space. 

Charge conjugation and Time reversal symmetries 



The Hamiltonian (S16l is symmetric under charge con- 



jugation (C) (particle-hole) transformation: 
U c c| CT U c = rna a 



(S25) 



where r\i is +1 on sublattice A and —1 on sublattice B. 
The unitary operator Uc can be explicitly written as 

U c = Y[e™ s *e^ G * (S26) 



where S°" are the spin operators acting on the singly occu- 
pied states and Gf are the pseudo-spin operators acting 
on the empty and doubly occupied states, and defined as 



G 



c t c l 



(S27) 
(S28) 



Every term in the Hamiltonian, Hq and T s , is C invariant: 



U C HU C = H 
U C T S U C = T s 



(S29) 
(S30) 



It then follows from equations ( S24 ) and ( S22 ) that every 
term of H k is C invariant, for all k. 



The time reversal operator is 



U T 4t U T = ia lo' C ™' 



U T = Y[e^JC 



(S31) 
(S32) 



where K, is the complex conjugation operator. Any state 
|hf) in the singly occupied subspace satisfies the condition 
G?\hf) = 0. It then follows that 



U C U T \M) = /C|hf) 



(S33) 



Time reversal symmetry of 

We show the TR symmetry of H^ k > by explicitly show- 
ing the equality of the matrix elements of and 
in a real basis. Specifically, we can choose 
the simultaneous eigenstates of Sf , 



S?\{<n}) = ^IW) 



(S34) 



We can always choose Sf to be real and hence we have 
K-\{ a i}) = K "*})- It then follows that, 



(W|C/- T ff<*>4|foy) = {{oiWU^UlUcH^UcUcUT^i} 1 ) 

= ({a^l/C^W/CK^}') = ({^}l^ (fc) l{^}'> (S35) 

I 

Thus the effective spin Hamiltonian is TR symmetric. Spin-Spin Correlation function 

This implies that it does not contain any odd-spin terms. 

We now outline the computation of the spin-spin cor- 
relation function. We write the Hamiltonian as 

H = H +H P (S36) 



■7i 



5f-Sf)+ 72 5fSf] 



(S38) 



J = 




*i = 



2i' 2 


3i' 2 






;^ 3 


t 2 - <' 4 


2C/ 3 



(S39) 
(S40) 



We take the Hamiltonian H (the Kitaev model) as the 
unperturbed Hamiltonian and Hp as a perturbation. We 
want to compute the correlation function 

ff(r,«) = (T(sS(t)5^(0))) (S41) 

Where r = r\&\ + r2e2, ei and e 2 are basis vectors as 



To leading order, this is the spin-spin correlation func- 
tion of the Kitaev model. The Kitaev model has a 6- 
spin conserved operator associated with every plaquette, 
W p which can take values ±1 and can be interpreted as 
a Z 2 flux. [7] The ground state is in the flux- free sector 
(W p = 1 Wp). The spin operators at site r create a pair 
of flux tubes in two of the plaquettes that the site be- 
longs to. Since the time evolution does not change the 



shown in Fig.(S5| and I, m are the sub-lattice indices. 



flux configuration, the spin-spin correlation ( S41 ) is zero 
except when r and are nearest neighbors. [8] 



The second-order perturbation term is 



,(2) 



H) 2 



dr\ I dT 2 (T 



(S42) 



The time evolution is governed by Hq. This term will be 
non-zero only if there are terms in Hp such that the prod- 
uct of the four operators in ( S42 ) do not change the flux 



configuration of the ground state. [5] We find that such 



terms do exist in H p 
tion {S^^ j^Sq q j^ 



We concentrate on correlation func- 
The following terms combine with 
A to produce flux-free configurations when acting 
on the ground state, 

"/2Sr 1 _ 1 ^ r2B S^ i _ lr2+1B ; l2Sl 1+lr2 _ 1A S^ 1+lr2A ; 

qX qV . qV qX 

H'~> ri ,r 2 ,A r 1 -l,r2 + l,B' < lJ r 1 ,r 2 ,A D r 1 -l,r 2 .B > 

These and the terms with (ri,f2) — > (0,0) which com- 
bine with Sq A give 36 possibly non-zero contributions 
to g( 2 \ 

The problem now is to compute the resulting 6-spin 
correlation functions in the Kitaev model. We do this in 
the fermionic representation of the spins in an enlarged 
Hilbert space, [7J in which the Majorana fermion opera- 
tors Cj and bf are defined as follows: 

of = {c i ,c j } = 25 ij (S44) 

{V? ) b$}=25 af) 6 ij , {ci,bf} = (S45) 

The physical subspace is defined by the constraint 

abp^ |^) phys = |^) phys (S46) 



In terms of these Majorana fermions, the leading order 
Hamiltonian is, 



Ho = J /J iCtCjibfb" 



(S47) 



This Hamiltonian describes Majorana fermions (cj), 
which we refer to as spinons, propagating in the back- 
ground of static Z 2 gauge fields (u^ a — ibfb") The 
correlation function in equation ( S42 1 thus factorizes into 



propagators of the Ci operators. Since the spin operators 
create two units of flux on adjoining plaquettes, the Ma- 
jorana fermion propagators are in the background of an 
even number of fluxes at a few points. 

We are interested in the asymptotic form of (r, t) 
when |r|, t — ¥ oo. Tikhonov et. al. [9] have shown that 
in this limit, the propagators are the same as those in 
the flux-free background. Their result can be physically 
understood by noting that the particles hopping far way 
from the flux pairs will not pick up any phases from them. 
Thus we can expect the long wavelength modes to be 
insensitive to a few localized flux pairs. 

To compute the asymptotic form or the propagators, 
we can derive the continuum theory of the low-energy 
modes in the flux-free background. The Hamiltonian 
(S47l, when ibfbj = 1, reduces to nearest-neighbor 



hopping on a honeycomb lattice just as in graphene. 
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Graphene has low-energy Dirac quasi-particles about two 
points, K and K', in the Brillouin zone. However, since 
the c fermions are Majorana fermions, the excitations ex- 
ist only over half the Brillouin zone. Thus the low-energy 
modes constitute a single Dirac quasi-particle. The con- 
tinuum theory is derived by introducing slowly varying 
fields ipi(r) such that 



Crl 



- 4K ' r W(r) 



(S48) 



Substituting equation ( S48 1 in equation ( S47 1 it can be 



seen that the low energy continuum theory is that of a 



single Dirac fermion. The propagator is defined as 

G lm (r,t) = (T(^(M)V4(0,0))) (S49) 

It can be computed to be 

1 1 



G hn = (r • r - Jtl)ir, 



47T ( r 2 _ J2{2)z 



(S50) 



where r = {r Xl Ty) are the Pauli matrices. We can then 
compute the correlation function in equation ( S42 1 to 
obtain the following expression: 



(SS(t)5&(0)) = (-0.56cos(2K • r) 7l 2 + 1.13 7 i 72 + 1.69 7 fe) detG + (0.56 7l 2 - 0.28 7 i 72 + 0M^ 2 e)Tr (Gt x G) 

+ (0.28 7l 2 + 0.07 7 | + 0.84 7l72 e + 0.63 7 2 e 2 - 0.28 7 i7 2 - 0.42 7 2 e)(TrG) 2 + 0.28 7l 2 Tr(r :r GT a; G) (S51) 



r 



where I represents the sublattice index (A or B), a rep- 
resents the three types of bonds x, y, z and e is the en- 
ergy density of the Kitaev model. We can obtain the 
(0)) and (S^(t)S^(O)) from the above corre- 
lation function by using the following property: a 2ir/3 
rotation about a sublattice A point takes x link to y link, 
y link to z and z link to x, in a cycle. The direction is 
reversed for sublattice B. Thus we can see that the ( S51 ) 



goes as r which shows that we have an algebraic spin 
liquid (ASL) for large U. 



Mean Field Theory of the effective spin model 



The effective spin Hamiltonian is given in equation 



(S36). To investigate the instability of the ASL to CSL, 



we perform a mean-field treatment of the above Hamil- 



tonian in the fermionic representation ( S44 1 . The decou- 



pling of the spinon and gauge field sectors is represented 
by 



-iciCj ibfbj 



P 



a/3 



The self-consistency equations are 



D 



a/3 _ 



ib?b 



C%j — (tCiCj 



(S52) 
(S53) 

(S54) 



We assume that the ground state is translationally in- 
variant, isotropic and denote 



a. 

B 



a/3 



n i,i+a. 



B 



= Wi Cj.iiei,; = m (S55) 

3U,i = fc; «^ ( S56 ) 



where a a represents the nearest neighbor vector on the 
a th link, a, (3 represent the links (x,y,z), I indicates 
the sublattice index (A or B) and e, represents the basis 



vectors of the underlying Bravais lattice (see Fig: S5 ) 





O 



FIG. S5: Basis vector used. Green (Blue) represents sublat- 
tice A (B). 



The mean field Hamiltonian at t' — 1 is, 
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Hmf = H MF + Hm F 



n MF 



u 



4 ( ( kl) ( k2) } I - iU ^ J UJ 



j e -iki 7 1 (e _ifcl +e ik2 ) 7i(e^ lfcl +l) 
U(k) = e I 7i(e- ifcl + e 4fe2 ) Je ik * 7i(e lfe2 + 1) 



L (e- lfel +l) 7l (e lfc2 + l) 



e ifc3 -e _ifel 
V 0i9 ,i = Mi72 | -e" lfc3 e lfc2 

— e — ** 3 —e~ ikl 
Vais.a = M272 I e lfe 3 -e"'* 2 



1 



keHBZ 



Vl(k) = 2«6i72 2J sm ( k ' e «) 

a 

u 2 (k) = -2i6 272 E sin ( k ' e «) 



J 



v <v _t w ™i( k ) * u ( k ) i f c ki 

4 . 4rL_ 1 kl k2 j I -i«*(k) zw 2 (k) J I c k2 



(S57) 
(S58) 

(S59) 

(S60) 

(S61) 
(S62) 

(S63) 
(S64) 
(S65) 



where h = k • e,- 




FIG. S6: Spinon dispersion relation at f7 = 2. 



The nearest-neighbor term in Hamiltonian ( S36 ) re- 



sults in closing the gap in the spinon sector at the Dirac 



points (see Fig: S6), whereas the next-nearest term col- 
lapses the gap at (0,0) and (it, it). For large values of 
U the nearest-neighbor term dominates; for smaller val- 
ues of U, on the other hand, the next-nearest neighbor 




FIG. S7: Spinon gap, E g as a function of U. 



term comes into play and the Dirac points shift to (0, 0) 
and (tt,tt). Numerically we find that the spinon sector 
is gap-less for U > 1.6 (see Fig S7 1 . We have checked 
that this remains true for 1 < t' < 0.5. VCA indicates a 
Mott transition at U — 2.4. This shows the absence of 
the CSL phase in the presence of higher order perturba- 
tive terms and indicates the ASL phase continues till the 
Mott transition. 
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